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ABSTRACT 



We study the dynamics of a twisted tilted disc under the influence of 
an external radiation field. Assuming the effect of absorption and reemis- 
sion/scattering is that a pressure is applied to the disc surface where the 



> 

cn 

00 

O^l ' local optical depth is of order unity, we determine the response of the vertical 



structure and the influence it has on the possibility of instability to warping. 

We derive a pair of equations describing the evolution of a small tilt as 
a function of radius in the small amplitude regime that applies to both the 
, diffusive and bending wave regimes. We also study the non linear vertical 

■ response of the disc numerically using an analogous one dimensional slab 

model. For global warps, we find that in order for the disc vertical structure 
to respond as a quasi uniform shift or tilt, as has been assumed in previous 
work, the product of the ratio of the external radiation momentum flux to the 
local disc mid plane pressure, where it is absorbed, with the disc aspect ratio 
should be significantly less than unity. Namely, this quantity should be of the 
order of or smaller than the ratio of the disc gas density corresponding to the 
layer intercepting radiation to the mid plane density, A <C 1. 

When this condition is not satisfied the disc surface tends to adjust so 
that the local normal becomes perpendicular to the radiation propagation 
direction. In this case dynamical quantities determined by the disc twist and 
warp tend to oscillate with a large characteristic period ~ X^^Tk, where 
Tk is some 'typical' orbital period of a gas element in the disc. The possibility 
of warping instability then becomes significantly reduced. 
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In addition, when the vortical response is non uniform, the possible pro- 
duction of shocks may lead to an important dissipation mechanism. 

Key words: accretion; accretion discs; binaries: close; galaxies: nuclei, x- 
rays: binaries, stars; hydrodynamics 

1 INTRODUCTION 

A significant number of X-ray binaries exhibit long-term periodicities on time-scales 
of 10-100 d. Examples are Her X-1, SS 433 and LMC X-4 see e.g. Clarkson et al. 2003. 
Precession and warping of a tilted accretion disc has been proposed as an explanation (Katz 
1973; Petterson 1975) which has subsequently been found to have observational support 
(Clarkson ct al. 2003). The effect of radiation pressure on a twisted tilted accretion disc was 
first considered by Petterson 1977a,b who noted that when radiation from the central source 
is absorbed at the disc surface and re-emitted, a potentially important torque will result. 
Iping & Petterson 1990 later suggested that such torques determined the shape of the disc 
and its precession rate. Pringle 1996 subsequently showed that an initially axisymmetric 
thin disc could be unstable to warping as a result of interaction with an external radiation 
field (see e.g. Maloney et a. 1998; Ogilvie & Dubus 2001 and Foulkes et al. 2006 for later 
developments) . 

The analysis considered the disc to behave as a collection of rings, which could interact by 
transferring angular momentum through the action of viscosity, but otherwise behaved as if 
they were rigid. In particular, the vertical displacement or tilt was assumed to be essentially 
uniform and independent of the vertical coordinate. Thus the pressure applied at the surface 
at optical depth unity is assumed to be effectively communicated through the disc vertical 
structure so as to give a near uniform response. In this paper we extend the theoretical 
treatment of the interaction of a disc with an external radiation field to take account of 
possible significant departures of the tilt or displacement from uniformity in the vertical 
direction. One of our objectives is to determine the conditions under which the assumption 
of a uniform response is valid, and then to estimate some of the consequences when they are 
not satisfied, including potential additional dissipation resulting from non linear effects such 
as the production of shock waves. The latter is done using a one dimensional slab analogue 
model with the required high resolution in the vertical direction. 

Although we focus on the effects of a surface pressure induced by an external radiation 
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field, very similar considerations are likely to apply when warps are induced by a surface 
pressure resulting from interaction with an external wind (e.g. Quillen, 2001) or surface 
forces produced by the interaction of the disc with an external magnetic field originating in 
the central star ( e.g. Pfeiffer & Lai, 2004). 

We begin by considering the relevant issues using simple physical arguments. Following 
Pringle 1996 we consider a thin disc immersed in an external radiation field with mid plane 
initially coinciding with a Cartesian (x, y) plane. It is supposed that the upper surface is 
parallel to the external rays so that there is initially no interaction with the radiation. The 
upper surface is then given a vertical elevation /i(r, 0), where we now use polar coordinates 
(r, 0). As a result, a disc element with surface area dA absorbs momentum from the radiation 
field at a rate 

Here Fq is the momentum flux at radius r associated with the external radiation field. The 
factor in brackets gives the angle through which the local normal is rotated as a result of 
the perturbation (see Appendix |A] for more details). Assuming the absorbed momentum is 
reradiated isotropically above the disc by the surface layers at an optical depth of unity, 
there will be an applied pressure there of magnitude 

This external pressure, when applied to a complete elementary ring, produces a net torque 
per unit length of magnitude 

where the complex vector 1 has components equal to (z, 1,0) in the Cartesian coordinate 
system. Here, in performing the azimuthal integration, we take into account that the az- 
imuthal dependence of h is through a factor of the form exp(— i0) and work with the radial 
amplitude from now on. 

To find the consequent evolution of the disc, one requires the component of the angular 
momentum per unit length perpendicular to its unperturbed direction. This is given by 

dj/dr = 27Tj:r^Q{h)i\, (4) 

with Q and S being the near Keplerian local disc angular velocity and the disc surface 
density, respectively, (h) = J dCph/T, where C is a vertical coordinate and p is the gas 
density. 
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For an elementary ring the condition that the rate of change of angular momentum equal 
the applied torque gives a tilt evolution equation of the form 



Assuming the vertical response is uniform, so that we can set h = (h) = rW, we obtain a 
description of warp evolution equivalent to Pringle 1996 when effects due to viscosity and 
bending wave propagation are neglected (see sections f5.ll and 1631 below). This description 
indicates the possibility of instability to radiation pressure warping. 

However, here we stress the fact that the averaged elevation (h), enclosed in angled 
brackets, applies to the bulk of the inertia of the disc and can therefore be shown to be close 
to the elevation of the mid plane. On the other hand h as used in the torque formula applies 
to the disc surface elevation. An important aspect of this paper is to distinguish these two 
elevations and investigate under what conditions they can be taken to be equal. This would 
be possible if the vertical structure responds as a rigid body to the external pressure forcing 
and we find the conditions for this to occur. The general requirement is found to be that 
the density at the surface of the disc where the pressure is applied should not be too small. 

It is possible to estimate in a simple way when the vertical displacement response of the 
disc to the external pressure forcing becomes non uniform. Let us suppose that the external 
pressure V is applied at the surface layer where the density p = p^. The induced vertical 
displacement /i, now should be considered to be also a function of the vertical coordinate C,. 

Assuming a linear response, which should be appropriate for sufficiently small elevations 
and, for simplicity an isothermal equation of state with sound speed c^, it can be easily shown 
that in the upper layers of the disc where vertical motions dominate the Lagrangian pressure 
perturbation has the form determined by presence of h, AP = — p^.c^||, see equation (HOj) 
below. Equating this to the external pressure we have 



If the vertical extent of the disc is (o, and h varies on a radial scale comparable to r, the 
characteristic change in h induced over the vertical thickness is easily estimated to be 



where we use the approximate relation = fi^o- 

From the condition Ah ~ h, using ([7]) and estimating |^ ~ /i/r, it is clear that whether 




(5) 




(7) 
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h is uniform or not is governed by the the parameter F^r / {pt,VL^C,Q){CQ/rY = A~^e5^, with 
A = p*/pc and 5 = (o/r so defining e = For /{pc^I'^(q). 

For a uniform response, one requires that <^ A. This is equivalent to the requirement 
that the product of the ratio of the external radiation momentum fiux to the local disc 
pressure with the disc aspect ratio should be significantly less than unity. We recall here 
that because of the geometrical configuration, the momentum fiux locally reradiated by the 
disc is in general much less than the external momentum fiux at that location. 

In this paper we extend the treatment of disc warping induced by the action of an external 
pressure to include the effects of the response of the vertical structure particularly when this 
is non uniform as is the case when the above condition is not satisfied. In that case we find 
that the disc dynamics enters a different regime. This is such that the exposed surface tends 
to align so as to reduce the momentum absorption and hence the applied pressure. In the 
extreme limit of this regime the upper surface acts as if it is in contact with a rigid wall 
with the tendency to radiation warping instability tending to vanish. In this limit quantities 
characterising the disc twist and warp tend to oscillate at typical frequency ~ XQ. 

To illustrate these effects we derive a description of the one dimensional evolution of 
the disc inclination in radius and time that incorporates the effects of the vertical structure 
response, radiation torques and which applies both to the high viscosity regime, when the 
evolution is diffusive, and to the low viscosity regime when the evolution is wavelike. In all 
cases, the efficacy of surface radiation pressure driven instabilities is found to be reduced 
once the model parameters are such that the vertical response is significantly non uniform. 

We go on to estimate conditions for the response to be nonlinear and investigate the 
development of shock waves in the response using a one dimensional slab analogue model 
which has the same characteristic behaviour of linear perturbations as the full disc model. 
We find that such shocks potentially provide an important dissipation mechanism. 

The plan of this paper is as follows. In section [2] we describe some aspects of the thin disc 
model used. In section [3] we go on to introduce the twisted coordinate system used together 
with the notation convention, giving the basic equations in section |H By integrating over 
the vertical direction we use these to derive a single equation governing the dynamics of a 
twisted disc in section |5l When the disc behaves like a set of rigid rings for which the vertical 
response is uniform, this equation can be used to give a complete specification of the warp 
evolution. We note that this provides an extension of previous formulations to be able to 
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consider the case when warps propagate as waves rather than diffuse radially (e.g. Nelson 
& Papaloizou 1999). 

In this Paper we use the twisted coordinate system formalism first introduced by Pet- 
terson 1977a, 1978 where the dynamical equations take the most simple form. When this 
formalism is adopted and an accurate description of all components of the equations of mo- 
tion is needed as in the problem on hand we show that, in general, there is an ambiguity 
in choice of twisted coordinates with a set of these describing the same physical situation. 
As discussed in section |5] a choice of a most appropriate twisted coordinate system can be 
motivated by the condition that perturbations of all dynamical quantities determined by the 
disc twist and warp are small. The transformation law between different twisted coordinates 
corresponding to the same physical situation is derived in appendix C for an inviscid disc. 

In order to obtain a complete description of the warp evolution when the vertical response 
is non uniform, we begin by obtaining a complete solution of the vertical problem for a 
polytropic model in section 16.11 This is used to obtain a pair of equations governing the 
radial warp evolution. We also indicate how the results can be extended to apply to more 
general models. We confirm the condition for the disc response to be like that of a series of 
rigid rings as A ^ e(^o/''^)^- 

We go on to perform a linear stability analysis of the radial evolution equations adopting 
a WKB approach in section [7] obtaining a maximum potential growth/decay rate in section 
17.31 In section [8] we discuss and confirm the analogy between the response of the disc and the 
linear and non linear dynamics of a vertically stratified one dimensional slab. We consider 
the development of shock waves and the formation of a rarefied hot atmosphere in section 
18.31 giving a crude estimate of the warp dissipation rate in section 18. 4[ Finally in section [9] 
we discuss and summarize our results. 

2 A THIN DISC MODEL 

To begin we briefly summarize some of the properties of an unperturbed axisymmetric thin 
accretion disc model which are used in the discussion below. 

However, the physical basis of our results does not depend on the detailed properties of 
this model which has been used for the purpose of making specific calculations. In order to 
make our treatment as simple as possible we use a highly simplified version of the Shakura 
& Sunyaev 1973 a-disc model. We assume a polytropic equation of state for the disc gas 
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where P and p are the gas pressure and density, r is the radial coordinate and 7 is the 
assumed constant specific heat ratio. It is assumed that the "polytropic constant" K is, in 
general, a function of the radial coordinate, r, but does not depend on the local vertical 
coordinate (, see equations (fT4l) . ( fTSl) below for a definition of the coordinate system used 
in our study. The equilibrium variation of density and pressure with height follows from the 
equation of state ([H]) together with the equation of hydrostatic equilibrium 



Here Cg is the adiabatic sound speed, Vt = ^ (GM) / is the Keplerian angular velocity and 
M is the mass of a central object. 
Integration of equation ([9]) gives 



P = Pc(r) 




1-7 





(10) 



VCo(r) 

where pc and Pc are the mid plane values of the density and pressure, respectively and Co is 
the disc semi-thickness. These quantities are related to each other through 



As shown below, many of our equations take an especially simple form for the case 7 = 5/3, 
accordingly we adopt this value below. 

It is assumed below that the standard Navier-Stokes equations are valid as the dynamical 
equations of our problem with the viscosity law set up according to the standard prescription 
(Shakura & Sunyaev 1973 ) 

z/ = an-^P, (12) 

where p is the kinematic viscosity. For simplicity we assume below that a is a small constant 
parameter and consider only leading terms in corresponding expansion series. 

The radiation intercepted by the disc is considered to result in an external pressure 
applied to a surface near to the free surface of the disc. This surface is formally defined 
through the condition p = p*(r = 1) where r is the optical depth. It is assumed hereafter 
that the ratio of the density at this surface to the midplane density p*/pc ^ 1- Thus, we 
neglect the effect of heating of the upper disc layers by incident radiation. Also, for simplicity 
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it is assumed that the radiation does not affect the structure of the unperturbed disc. This 
can be achieved by setting the disc semi-thickness to be proportional to the radial distance 

Co = Sr, (13) 

where 5 -C 1 is the disc opening angle which is constant. 



3 COORDINATE SYSTEM AND NOTATION CONVENTION 



The dynamical equations governing a twisted disc take their simplest form in the so-called 
twisted coordinate system first introduced by Petterson 1977a. Since this system has been 
discussed extensively elsewhere (for details see e.g. Petterson 1977a; 1978), in this paper we 
give only some basic definitions and relations important for the future discussion. 

In order to define the twisted coordinates let us consider a fixed Cartesian coordinate 
system {xi, X2, x^) with origin at the central star and (xi, ^2) plane coinciding with the mid- 
plane of the unperturbed disc. We then perform two elementary rotations of the coordinate 
axes parametrised by the two Euler angles 7 and f3. The first rotation is of the (xi,X2) 
axes through the angle 7 keeping the X3 axis fixed. As a result of this rotation the Xi and 
X2 axes change their directions. The second rotation is of the (2:2, ^3) axes in the plane 
perpendicular to the new fixed direction of the Xi axis and is through the angle The 
cylindrical coordinates (r, ip, Q defined with respect to the rotated axes now determine the 
twisted coordinate system. 

Formally, in the linear approximation that the inclination /3 is small, the transformation 
to the twisted coordinate system {r,ipX) from the initial Cartesian coordinates {xi,X2,X3) 
is given by 



(14) 



where B is the rotation matrix 

/ cos 7 sin 7 \ 

B{I3,^)— — sin7 COS7 (3 , (15) 

\ /3 sin 7 —(5 cos 7 1 / 

and f3{t,r) and '~f{t,r) are, in general, functions of r and the time t. 

Instead of the original Euler angles it is useful to introduce the combinations ipi = (3 cos 7, 

ip2 — Psinj. Further, we define the angle (p — ip + j which is used to replace the angle 



/ rcos-^X 






r sin'^ 


= ^7) 


X2 


V c J 




U3/ 
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if). We are thus able to consider the case f3{r,t) = for which the transformation (fT5|) 
is degenerate. The rotation matrix giving the transformation from the twisted coordinates 
back to the Cartesian coordinates is given by transpose of B, . The cyhndrical coordinates 
{tc, 0C5 -2c) associated with the original Cartesian system are defined through: Xi = r^ cos 0c, 
X2 = sin 0c and = Zc- From equations ([T^ and f[T5l) we obtain 



sin0c = sin0 — — cos0cos(0 — 7), (16) 
r 

rc = r- PC sin(0 - 7) , (17) 
and 

Zc = C + Prsm{4)--f). (18) 

From equations ( JT6l) and (fTTll it follows that when ^ = we have (f)c = 4> and Tc = r. In 
addition, the difference 0c~0 is also first order in the small angle (3. Thus we can set (pc = 4> 
in equations for quantities that are first order in (3. 

All vectors and tensors appearing in our calculations are projected onto a local orthonor- 
mal basis e^, e^, and (Petterson 1977a, 1978) where the bar is used to denote a vector 
quantity from now on. The explicit expressions for the above basis vectors are not important 
for our purposes and can be found elsewhere, see e.g. Petterson 1977a, 1978, Demianski & 
Ivanov 1997 hereafter DI. We note, however, that the basis can be uniquely defined by the 
requirement that e,^ is the coordinate vector such that ■ V = ^ for C ^ 0. All other 
vectors can be found from the orthonormality conditions together with the requirement that 
the coordinate axes form a right handed system. 



3.1 State variables 

Since the coordinate lines of the twisted coordinate system are not orthogonal and both 
the coordinate lines and directions of the basis vectors depend on time and on the radial 
coordinate, differentiation of quantities is best performed with help of appropriate connection 
coefficients. Their explicit form can be found in Petterson 1977a, 1978 and DI. Also, we 
note that the projections of the velocity vector v onto the basis vectors are not in general 
proportional to the time derivatives of the appropriate coordinates of a particular fluid 
element, see Petterson 1978. Accordingly, we give relations between these components and 
the time derivative of the appropriate gas element coordinates where this is explicitly needed. 
Writing them in the twisted coordinate system, in the approximation that the inclination 
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angle is small enough that a linearization procedure in which quadratic and higher powers 
of P can be neglected, the Navier-Stokes equations and the continuity equation then lead to: 

1) a set of equations for unperturbed quantities Qo formally coinciding with the equations 
describing them for flat disc accretion. 

2) A set of equations for quantities Qi, which are non axisymmetric perturbations to Qq 
induced by deformation of the disc shape. The dependence of the latter quantities on (p is 
harmonic: Qi = Qi'^ cos (j) + Q^^ sin (f) Q 

It is convenient to represent these perturbations as well as the angles jS and 7 using the 
complex quantities, Q and W, which are defined through 

Q = QS'^ + ^Q!'\ W = f3e'^ = ^1 + t^2. (19) 

Here we note that we adopt the convention that quantities in bold represent complex 
perturbations. This convention applies also to the components of a vector with vectors 
themselves denoted with an over bar. When real as in the perturbed equations of motion 
(I23l)-(l29l) given below, these are not in bold. Note also that for convenience we use pi and 

to denote the complex density perturbation and vertical component of the Lagrangian 
displacement below, and these are not in bold. 



4 BASIC EQUATIONS 

As we mentioned above we write each of the variables entering the continuity and the Navier- 
Stokes equations as the sum of an unperturbed part, formally coinciding with a solution of 
the equations for axisymmetric thin disc accretion, and a perturbation. Accordingly, we have 

v' = vi, + vi, p = po + pi, (20) 

for the components of velocity and the gas density respectively. But as the only non zero 
component of the unperturbed velocity is in the azimuthal direction and corresponds to 
Keplerian rotation which we deal with explicitly, we shall drop the subscript 1 from the 
velocity perturbation. 

The unperturbed gas density po is given by equation ffTOj) . We assume that only circular 
Keplerian motion of the gas is present in the unperturbed state and accordingly have Vq = 

^ Strictly speaking, in order to satisfy boundary conditions at tiie disc surface we can assume a harmonic dependence of some 
of these quantities only over a limited domain of the angle <f>, see Appendix |A] for details. However, as follows from the results 
provided there, we can formally proceed assuming such quantities depend harmonically on (j) over the appropriate domain and 
then extend to the whole <f> domain using the rule discussed in Appendix IaI 
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rQ = y^^- Thus we neglect the radial component of the unperturbed velocity Vq present in 
the thin disc mainly due to action of viscosity. This can be done in the limit of a sufficiently 
small value of a < 1. 

Apart from two significant modifications, our equations for the perturbed quantities 
follow from the corresponding equations given in DI provided that terms involving quadratic 
and higher order corrections in the assumed small parameter a are discarded and that the 
gravitomagnetic force in the ( component as well as relativistic post-Newtonian corrections 
to the r and 0-components of the Navier-Stokes equations are neglected. 

The most important modification stems from the fact that in the present study an accu- 
rate description of vertical motions and displacements is needed and therefore, we take into 
account all contributions to the vertical velocity, v'^. These consist of the contribution deter- 
mined by time-dependent evolution of the basis vector and the contribution determined 
by time dependent vertical displacements of gas elements in the disc. Accordingly, we have 

v'^ = rU + v'^. (21) 

where lA = (5 sin ip — ^'j cos = sin (/) — \&2 cos (j) and the dot stands for the time derivative. 
The relation of f to the corresponding component of the displacement vector ^, is given 
by Petterson 1978 

v^ = i' + n^^'- (22) 

Only the first term in (PT]) . which arises from the time dependence of the changing coordinate 
system, is normally taken into account in studies of twisted discs that use the formalism 
of the twisted coordinate system, see Hatchett at al 1981 and DI. However, the role of the 
second term, which describes the disc motion relative to the moving coordinates, is essential 
for our purposes. 

We also take into account the time derivative of the perturbed density pi in the continu- 
ity equation and the time derivative of in the (^-component of the Navier-Stokes equation. 
These terms are proportional to a small parameter l/{Qtev) where tev is a characteristic 
evolution time of the twisted disc. That is for any perturbation quantity Q, we can write 
IQ/QI = 0{l/{Qtev))- However, even though they are small, in a formal sense, when com- 
pared with leading terms in a series expansion, when retained, they allow us to regularise 
some otherwise formally singular expressions, that arise when dealing with the particular 
problem of calculating the vertical displacements in the disc, see e.g. equation (!50|) below. 
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The perturbed Navier-Stokes equations can be divided into 'horizontal' part incorporat- 
ing the (r) and (0)-components and 'vertical' part corresponding to the {() component. The 
horizontal part follows from the equations provided in DI when the relativistic corrections 
are neglected and the limit of small viscosity is adopted. The (r)-component can be written 

as 

Po{v^ + n{^v^- - 2v^)) = -po^^VCW - (23) 
where W = f3' sin(-?/') — (3'j' cos ip = sin — \l/2 cos and 

t< = -ri-^v^ (24) 

is the (r, C)-component of the viscosity tensor. In the expression for W, a prime denotes 
differentiation with respect to r. 

Differentiating the (0)-component of perturbed Navier-Stokes equations with respect to 
(0) we obtain 

where 

t^' = -V§^vt (26) 

The set of equations fl^Tl) and fl2^ can be used in order to express the perturbed velocity 
components in terms of the quantity W. As seen from equations (12 ip and (1241) . to the 
leading order in small parameters a ( or equivalently 77) and l/{Qtev), the perturbed velocity 
components enter in both equations in the same combination 

|.-2.*. (27) 

Therefore, these equations are degenerate to leading order and accordingly the next order 
corrections have to be retained (Papaloizou & Pringle 1983, hereafter PP). This leads to an 
inverse dependence of the 'horizontal' part of the perturbed velocity on the small parameters 
such that both and v"^ are oc min(a^-'^, 1/ (Qtev)) (see PP for details). 
The (C)-component of the Navier-Stokes equation is written in the form 

Po + n^v'^ + 2nr^U^ = -n\p, - ^A, (28) 

where pi and Pi are the perturbations of density and pressure, respectively. These are related 
through Pi = c^pi with the sound speed given by equation ( |TO|l . 

The continuity equation is formally equivalent to perturbed continuity equation written 
for the case of a flat disc in cylindrical coordinate system 
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(29) 



We recall that the relation between the velocity component and the corresponding com- 
ponent of the Lagrangian displacement vector ^ is given by equation fl22|) . 

Apart from the terms proportional to W and U equations ( |23ll29l) are formally equivalent 
to the corresponding equations that would be written down directly in a fixed cylindrical 
coordinate system (r, 0, ^) (PP). 

The term proportional to W in equation (1231) is essentially a projection of the pressure 
gradient. Due to the presence of perturbations with odd symmetry with respect to ( in 
twisted discs, surfaces of constant pressure do not in general coincide with the orbital planes 
of gas elements. Therefore, there is is a projection of the pressure gradient onto such an 
orbital plane which is described by this term, see Fig. 4 of Ivanov & Illarionov 1997 for 
graphic representation of this effect. 

The term proportional to in equation (128!) is made up from two contributions. The 

first is related to the non-holonomic character of our basis vectors. As we have mentioned 
above in a non-stationary situation, this leads to a non-zero vertical velocity even for a 
gas element being at rest with respect to the twisted coordinate system, see equation ( l2T|l 
above. This fact was noted for the first time by Hatchett, Begelman & Sarazin 1981. The 
second contribution is determined by non-inertial effects. These are taken into account by 
an appropriate connection coefficient in the formalism developed by Petterson 1978. Both 
contributions have exactly the same form. This accounts for the factor 2 in the expression 
2^W in equation (|28D . 

The set of equations fl23tl29p can be brought into a simpler form with help of the complex 
notation introduced through equation (fT9!) for the density perturbation and the perturba- 
tions to the velocity components. Also, we assume hereafter that all perturbed quantities 
depend on (p and time only through an exponential factor. Thus 



where the frequency oo is, in general, a complex constant. 

Subtracting equation ( l25i) from equation ( l23l) and introducing the complex notation it 
is easy to see that the result contains the components of the velocity perturbation only in 
the combination v+ = v*" + 2^*^. The resulting equation can thus be represented as 



Q oc e 



(30) 



dW 



(31) 



dr ' 
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where we use equations (12^ and (126!) and set p = po from now on. The quantities v*" and 
v*^ can be separately expressed in terms of v+ taking into account the fact that in the 
leading approximation in the small parameters we can equate expression ( |27ll to zerj^ and, 
accordingly, obtain 

v*^ - 2zv'^ = 0, = -v+, v"^ = --V+. (32) 

2 4 ^ ^ 



Now we can express v'' and v"^ in terms of W with help of equations ( l3T|) and ( l32l) . To do 
this we remark that the ( dependence is dealt with by noting that r] oc P and the solution is 
such that both v'' and v*^ are oc (. We then use the resulting expressions to eliminate these 
velocity components in the continuity equation fl29l) so obtaining 

{u + Q)p, = -I Lac + ^(pvO j , (33) 



where 

7-1/2 Qr \ Co dr 



A = -tttP V — r 7- , (34) 



and 

il) = uj + iaQ. (35) 

Finally, equations fl28|) and fl22|) can be brought in the form 

d 

- ip{uj + Q)^r^ = 2ipQujrW - Q^Cpi - g^clpi, (36) 
and 

= -i{uj + fi)e^. (37) 



5 DERIVATION OF A SINGLE EQUATION DESCRIBING THE 
DYNAMICS OF A TWISTED DISC 

At first let us point out that the hydrodynamical equations as written in a twisted coordi- 
nate system cannot in principle provide a single dynamical equation for W without some 
additional considerations. Indeed, since the Euler angles (3 and 7 are considered as dynami- 
cal variables in this formalism, the total number of dynamical variables is always larger by 
two than the number of equations actually needed to determine the dynamics of the system. 

^ Setting II27I I to zero physically reflects the fact that trajectories of gas elements in the disc may be approximated as Keplerian 
ellipses in the leading order. Next order correction determine parameters of the these ellipses and effect of slow precession of 
their main axes on time scale of order of tev, see DI and Ivanov & lUarionov 1997 for more details on physical picture of the 
horizontal motions in the disc. Formally, parameters and precession rate of these ellipses can be obtained from equation II31I I. 
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For example, in the case of a baratropic equation of state P = P{p), there are six 
dynamical variables used in the description of the problem, these being the three components 
of velocity, and the density together with the two Euler angles. However the number of 
dynamical equations available is four, these being the three components of the Navier- 
Stokes equations, and the continuity equation. Since the number of equations is smaller 
than the number of variables, there are in principal, different sets of variables v{t, r), pi(t, r) 
and W(t,r) describing the same dynamical system and connected to each other by certain 
transformation laws. These transformation laws can be easily found from the condition 
that the velocity and density fields measured in some inertial coordinate system remain 
unchanged when the set of dynamical variables associated with the twisted coordinate system 
is transformed. Explicit expressions for these transformation laws are given in Appendix O 
It is interesting to note that the dynamical variables in the inertial coordinate system play 
the role of gauge independent quantities and the laws of transformation between different 
sets of variables in the twisted coordinate system may be regarded as gauge transformations. 
This situation has many analogies in different physical systems, say in systems governed by 
General Relativity, such as e.g. , the theory of small cosmological perturbations, see e.g. 
Landau & Lifshitz 1975. 

In order to begin the process of obtaining a single equation describing the dynamics 
of a twisted disc referred hereafter as a 'dynamical equation' we substitute the density 
perturbation given by equation ( l33i) into the 'vertical' equation ( 136|) thus having 



where we use the hydrostatic balance equation ([9]) and approximate Q + u ^ Q in the 
first term on the right hand side. Equation ( l38l) plays an important role in our analysis. A 
dynamical equation follows from ( l38i) provided that appropriate boundary conditions at the 
disc surface are specified. 

5.1 A simple approach to the derivation of a dynamical equation 

Before discussing an approach to the problem which fully takes into account the vertical 
structure of the disc, we would like to show how to obtain a governing dynamical equation 
from equation fl38|) in a straightforward way. Although this approach is incomplete it allows 
us to obtain a dynamical equation which is approximately correct provided that the ratio of 




d_ 



(c'poAC), 



(38) 
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the density where the external pressure is apphed to the central density, A = p*/pc, is not 
too small. 

To do this we integrate equation (138|) over ( taking into account the fact that the last 
term on the right hand side vanishes when C ±(o and can be discarded. We thus obtain 

(^c^p^v? + clpAC^ \^-% + 2QUJ pv^rfC = 

. GMQ'^ d d \ 
-2T.n^ujrW T^^\ rW , (39) 

2rV2 dr\ uj dr I' ^ ' 

where S = dC,po is the surface density and a typical disc height H is defined by condition 
= c^CCVo/S. 

It can be easily shown that the surface terms on the left hand side are, in the low 
frequency limit, proportional to the Lagrangian pressure perturbation, AP = [dP^/dQC,'^ + 
Pi. Introducing the complex notation and taking into account that in the low frequency 
limit, we can write from ( 1371) that ~ ifi^^v'', we obtain 

AP - + c^Pi - (c^P^v^ + c^pAc) , (40) 

where we make use of equation (|33l) and use the hydrostatic balance condition ([9]). 

The Lagrangian pressure perturbation must be equal to the radiation pressure at the 
surface corresponding to an optical thickness expected to be of order unity. As we discuss 
in Appendix [A] we can assume that the radiation pressure acts only on the upper surface of 
the disc ( i.e. where ( +Co)- We set, accordingly, the surface term corresponding to the 
lower surface of the disc to zero, and express the term corresponding to the upper surface 
in equation fl39|) with help of equation fHOl) through the radiation pressure term F_|_ given 
by equation ( lAllI) of Appendix [A] to obtain 

where we neglect the second term of the left hand side of ( l39l) proportional to uj in order to 
get (HTi) . This procedure is justified below, see Section 6.4 where we derive a similar equation 
in a more rigorous way. 

Also, let us stress that the (C)-component of the displacement vector is assumed to be 
evaluated at the surface of the unperturbed disc where the optical thickness r is of the order 
of unity, and let us recall that Fq is the radiation momentum flux per unit area perpendicular 
to the radial direction, see Pringle 1996. 
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The density p*(r ^ 1) where the external radiation pressure is apphed, is assumed to be 
always much smaller than the mid plane density pc- 

We comment that in the absence of radiation pressure, equation (l4Tll describes the linear 
dynamics of a free warped disc. When a = 0, this is wavelike, the waves being non dispersive 
and having local speed QH/2, see Nelson & Papaloizou 1999, Papaloizou & Lin 1995, Nelson 
& Papaloizou 2000. Furthermore we comment, without giving details, that equation (HTl) may 
also be derived following the approach discussed in those papers, so confirming the analysis 
based on the adoption of a twisted coordinate system presented here. 

The radiation pressure term involves the combination C.^'/r + zW. This can be regarded 
as being the angular displacement of the surface of the disc. It is composed of two parts. 
The first ^^/r represents the displacement of the surface of the disc relative to the midplane. 
The second term zW can be viewed as representing the displacement of the midplane and 
is described by the twisted coordinate system. When the disc behaves like a rigid ring, the 
contribution of is small. 



5.2 The small parameter A = p^,/pc 

The small parameter A = p*/pc should be compared with other small parameters of the 
problem in order to specify the different possible regimes of dynamical evolution of the disc. 
In particular, as we show below it is important to compare it with the local growth rate, 
in units of Q, associated with any possible instability of the disc that is driven by radiation 
pressure, Ur- 

As mentioned in the introduction and elaborated further below, we will see, when A > ujr, 
we can neglect C,'' in the expression .^^/r + iW entering in (HTl) . Physically this corresponds 
a situation when the surfaces of constant density in the disc are approximately parallel to 
each other and remain at an approximately constant distance apart. In that case, for the 
purpose of calculation of the radiation pressure term, the disc may be considered as being 
composed of rigid rings having different inclinations and precession angles and the analysis 
undertaken in previous studies (e.g. Pringle 1996) applies. Neglecting the displacement 
in dm), we see that this equation gives the evolution of W without the need for additional 
considerations. For the case a > 6 the frequency uj may be neglected in the expression 
u = u + iafl and equation is reduced to a form analogous to that obtained by Pringle 1996. 
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In this case Pringle 1996 indicated that the disc may be unstable in hnear theory provided 
that the radiation pressure term is sufficiently large. 

In the opposite limit 6 > a, equation ( 1411) describes wave-like propagation of twisted 
disturbances of the disc (see Nelson & Papaloizou 1999, Papaloizou & Lin 1995, Nelson & 
Papaloizou 2000) as indicated above, modified by the presence of the radiation pressure 
term. 



5.3 The polytropic model and the dynamical equation in a 'standard' form 

For the polytropic equation of state we can bring equation (14T!) to a further reduced form. 
In this case one can obtain from equation ffTOl) 

where T{x) is the gamma function. 

Before derivation of the general dynamical equation let us now temporary assume that 
the disc parameters are such that A > 6"^ /a and the term proportional to in fHTl) can be 
neglected and obtain an equation analogous to what has been obtained by Pringle (1996) 
for the polytropic model. In this case we can use equation (142|) and rewrite equation (14T!) in 
the form 

where and introduce the dimensionless radiation pressure parameter 

^" ^ Jc! 

defined as the ratio of the radiation momentum flux Fq to a characteristic gas pressure in 
the mid plane of the disc. In general, the condition Fq <^ 1 should be fulffiled for a simple 
description of the evolution of a twisted disc to be valid. Let us note that when 7 = 5/3 
the numerical factor in the second term on the left hand side of (l43l) is equal to ^ and the 
numerical factor in the last term is equal to ^. 

Equation fj43|) describes the 'standard' dynamics of a twisted disc immersed in a radiation 
field. 
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6 THE SELF-CONSISTENT APPROACH 

Now let us take into account the contribution of all the terms in equation (|4T1) . In order 
to do that we need to relate the term in this equation involving the {() component of 
the Lagrangian displacement vector to W. Also, the term 2fluj J^^ pv^d( was discarded 
without justification in order to go from equation (1391) to equation (14T]) . We provide a 
formal justification for this below. In order to find the vertical component of velocity, and 
accordingly, the vertical displacement we solve equation (138!) considering the right hand side 
of this equation to be a source term. From the appropriate solution of this equation, we show 
that the dynamical equation has indeed the form of equation (jUj) to leading order in our 
small parameters a and l/{Qtev) and we also specify the form of the vertical displacement 
at the surface where the external pressure applies. 

6.1 Solution of the vertical problem for the polytropic model 

In order to find a solution of ( |38l) it is convenient to separate the vertical velocity into 
two parts such that v*" = Vq + v^. 

Here the first part of the decomposition, Vq, is a particular solution of (1381) with boundary 
conditions at the disc surfaces which apply when they are free and which are regular in the 
limit of vanishing density. They are such that in that limit, the Lagrangian perturbation of 
the pressure given by equation (l40l) should approach zero at the disc surfaces ( = ±Co- 

The second part of the decomposition, v^, is a solution of the homogeneous part of 
equation (l38l) : 



Yc ("'^^"^O ^ ^^^""^^ ^ °- ^^^^ 

This solution will be used in order to satisfy the boundary conditions 

AP(p = p*, C ^ +Co) = F+(p = p„ C ^ +Co), AP(p = p*, C -Co) = 0, (46) 

where the Lagrangian perturbation of pressure is given by equation fHOj) . We comment that 
expressions such as C ^ ±Co are used below to apply to the surface where p = p*, the latter 
being the small density where the external pressure is applied, rather than at strictly zero 
density. 

Firstly we derive an expression for the quantity Vq. To do this we transform the right 

hand side of equation fl38l) with help of equations ffTOl) and ffTTj) to the form 

f) CM 
S ^ -2pQ'u;rW - pQ'AC' - ^(c^pAC) = -p^^oorW + (7 - 1)^B), (47) 
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where 

B = CoV;'^ (PcN) - ^x\rp:'^ (Co""^PcN) , (48) 

X = C/Co and N = r"'^^"^ [dW / dr) / Cj . Note that we do not assume that the disc opening angle 
is constant when deriving the relations P7|) and (H8l) and these are valid for any dependence 
of Co on r. 

The calculation of Vq is greatly facilitated by noting that the source term S/p has a 
quadratic dependence on x. Is is also easy to see by direct substitution that a solution for 
the vertical velocity Vg with the same quadratic dependence on x can be readily found in 
the form: 



bi + h2x\ (49) 



where expressions for the quantities bi and b2 directly follow from equation (155]) and the 
form of the source term ( H71) 

and 

3 (7-1) GM^^ d / \ 



4 (37 - 1) ri/2 

where we introduce the dimensionless frequency uj = u/Q and assume that it is small: 
|cl;| <^ 1. In this limit it is easy to see that the quantity b2 is ~ a; times smaller than bi. It 
is, therefore, neglected later on, and we assume that 

^ bi (52) 

in our future analysis. 



6.2 The solution of the homogeneous problem 

Now let us find the homogeneous component v^. It turns out that the solutions of equation 
(145|) take on a very simple form when 7 = 5/3. In that case the solutions can be expressed 
in terms of elementary functions. Therefore, we shall specialise to this case in the main 
text of the paper from now on. However, the final description of the system we obtain 
has wider applicability (see below). The form of the solutions of the homogeneous problem 
corresponding to general values of 7, which can be used to construct the full solution in that 
case is relegated to Appendix [Bl 
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When 7 = 5/3 it is convenient to introduce new variables 6 and Y according to the 
prescription 

C = Co sin e, = Y/ cos^ 6. (53) 

It is obvious from equation (|53|) that ^ = and 6 = 7i/2 correspond to the disc midplane 
( = and to the free surface C = Co; respectively. It follows from equation (ITOl) that 
the density is proportional to cos^ 6 and accordingly the variable Y is proportional to (- 
component of the gas momentum per unit volume. 

In terms of the variables defined through (l53l) . equation fj45l) takes the form 

-^Y + 2tane-^Y + 3(l + 2cD)Y = 0. (54) 
do do 

It can be verified by direct substitution that the two independent solutions of (15^ are 
Yi = (k+1) cos(/s:-1)^+(k-1) cos{k+1)9, Ya = (k+1) sin(/«-l)^+(K-l) sm{K+l)9, (55) 



where k = V4 + 6uj. Note that the solutions ( l55l) are valid for any value of u. For our 
purposes, however, only the case ^ 1 is significant and we accordingly set k = 2 + |a; 
and take into account only zeroth and first order terms when summing series in ascending 
powers of c^. 

Thus the general solution of equation (J38l) can be expressed as a linear combination of 
these independent solutions in the form 

v< = (C^Y^ + C2Y2)/ cos=^ 9 + (56) 

where is given by equation f j49l) with 7 = 5/3 and the arbitrary constants Ci and C2 can 
be chosen to satisfy the boundary conditions fH6l) . 



6.3 A specific gauge freedom associated with the twisted coordinate system 

Following the simple approach of section 15.11 a dynamical equation for the quantity W 
can be obtained from equation (!38l) by integration over ( provided that certain terms are 
discarded. However, in a more exact approach, it is easy to see that this equation cannot be 
used for determining the evolution of W without invoking some additional considerations. 
As we have mentioned above this situation comes about because there is some degree of 



^ Contrary to the inhomogeneous part Vp the next order terms in w in the homogeneous part diverge near the disc surface. 
They can be of order of or larger than the leading terms and must, therefore, be retained. 
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arbitrariness associated with the specification of the twisted coordinate system used in our 
study. 

The resulting gauge freedom may be used in order to put constraints on some dynamical 
variables by choosing the most appropriate gauge from physical point of view. In our case 
it seems reasonable to look for a gauge where the vertical component of velocity has a 
small absolute value. This could be specified by imposing the requirement that 

.;^(C = 0) = 0. (57) 

As we will see below the condition f lFTI) together with equation (1551) fully specify all dynamical 
variables and allow us to obtain a single equation determining the dynamical evolution of 



W 



.0 



6.4 Derivation of the dynamical equation 

In order to find the dynamical equation we calculate the coefficients Ci and C2 entering (!56ll 
with help of equations ( l40l) . ( H6l) and ( lAllI) . Then application of the gauge condition ( 1571) 



will give us the dynamical equation. 

We note that only the homogeneous part of the vertical velocity can lead to a non 
vanishing Lagrangian pressure perturbation in the limit of zero density when ( ^Co- Thus 
this component of the solution dominates near the surface. Using equations (HOl) . (H6ll and 
( lAllI) . in the limit of small surface density, p*, we get 

where we assume that all quantities are evaluated a.t p = p^. 

Differentiating expression (156|) , using (155|) remembering that C = Co sin ^ and considering 
the limit ( ±Co and accordingly 6 ±7r/2 we approximately have 

^AiC - ±Co) ^ l^^iT^^^C^ + 2C,). (59) 
o( Co cos*^ 2 

The second condition in (155]) tells that C2 = —\tiCjCi. Using this fact we obtain from 
equations (l52l) . (155|) . (l56l) and (l59l) in the same limit 

|:V«(C - +Co) ^ ^MiC - +Co) ^ -T^C^i, (60) 
dC dC Co cos^ e 

* Note that the condition Il57|l is not unique. Another reasonable condition would be eg. the requirement that J^^° <^CP?C ~ ^' 
We expect, however, that any reasonable condition used to fix the gauge would give the same dynamical equation to leading 
order in the small parameter l/{Qtev)- 
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and 

= 4 + ^ b, + Ci(4 - Ip^) = b, + C,(4 - (61) 

where we recall that A = p^:/pc = cos{d) such that the coordinate 6 corresponds to the 
density level where p = p^. 

Using equations ffTOj) and ffTTl) with 7 = 5/3 we can represent the factor in the front of 
the radial derivative on the right hand side of equation (!58|) in the form 

l^^r = ^nr. (62) 
3 cjp^ cos^ 

Finally, substituting equations (160|) and (16T|) in equation (158!) and taking into account (162|) 
we obtain an equation for the quantity Ci 

where we introduce an important parameter 
Co '~ Pc^Kl 

This parameter determines strength of dynamical effects caused by the radiation pressure 
acting on the disc. When e > 1 these effects may be significant. Physically, the inequality e > 
1 means that the ratio of the momentum flux carried by radiation, Fq, to the characteristic 
mid plane gas pressure in the disc is larger than the disc opening angle. 

Now let us consider the gauge condition fl571) . From this condition and the expression 
fl49|) it follows that 

C,Y,{e = 0) = -bi, (65) 

where we use the variable Yi defined in equations fl53|) and fl55|) . Taking into account that 
Yi (6* = 0) =4 and, accordingly bi + 4Ci = in equation (l63l) and using the explicit form 
for the quantity bi given in equation fl50l) with 7 = 5/3 we obtain 



i = -h = -rr^. (64) 



Equations (1631) and (1661) form a complete system for describing the disc dynamics. Substitut- 
ing Ci given by equation ( l66l) into equation ( l63l) we can obtain a single third order equation 
for the quantity W. However, we find it more convenient to retain the pair of equations 
obtained directly from (163!) and (166|) to describe the disc dynamics. 
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6.5 Description of the dynamical system 

In order to bring this pair of equations into a simpler form, we introduce a new dimensionless 
variable x = ^Ci. Using x to replace Ci in equations (!63|) and (!66|) we obtain 



3 

-X = 
2 

and 



(67) 



x-^W=<^AfC|^i^wU(^i^f5^AwV (68) 

24CoPc dr\ u dr J A'Z dr \ Cj dr j ' ^ ' 



where we have used the explicit expression for the quantity N defined through equation fl48|) 
and recall that Co = uo + iaVt = Q{uj + ia). 

It is easy to see that equation ( 14T|) follows from form equations ( l67|l and ( l68i) provided 
that the ((^)-component of the Lagrangian displacement vector in equation ( HTi) is replaced 

by 

e« = - — A-Vx. (69) 



In order to be able to neglect this term in equation fl67j) . we require that |x| ^ |AW|. 
When this is satisfied, x can be determined after neglecting the first term in this equation, 
and assuming that W varies on a length scale comparable to r, with the result that x ~ e5^. 
Thus a description assuming that the vertical response of the disc is rigid can be recovered 
only when the disc parameters are such that A ^ e5^. This is precisely the condition we 
obtained from a simple argument given in the introduction. 

In the opposite limit the contribution proportional to ^(A~^x) should not be neglected 
in equation ( !67l) when the efi^ects of the radiation pressure are significant. Then as A — >^ 0, 
this equation gives 

X = i-AW. (70) 

Using the above to eliminate x in equation fIBSl) we obtain an equation describing stable 
radially propagating disturbances which attain a frequency given by = -^A in the limit of 
zero radial wave number (see also the discussion below). 

We comment here that although we used a polytropic model to derive equations fl67|) 
and (l68l) and continue with the immediate discussion, the indicated equivalent description 
of the dynamical system, which uses S and H with the numerical factor Sit /A replaced by 
2S/(pcCo)? contains no explicit dependence on the local details of the model, other than 
through the product of the height at which the external pressure is applied and the value of 
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the density there, S and H. This shows that our formahsm has a more general apphcabihty 
and we have verified this to be the case by obtaining the equivalent representation, quite 
generally, using an alternative analysis that does not use twisted coordinate systems, see 
e.g. Nelson & Papaloizou 1999, Papaloizou & Lin 1995 Nelson & Papaloizou 2000. 



7 QUALITATIVE ANALYSIS OF THE DYNAMICAL EQUATION 
7.1 Local analysis 

In general, equations (1671) and (|68l) should be solved numerically, for a specified disc model. 
However, this is beyond the scope of the present paper. Here we follow Pringle 1996 and adopt 
the usual local approximation scheme in which it is assumed that the radial dependence of 
all variables is cx e^^^ and that the radial wave-number k is such that \kr\ ^ 1. Note however, 
that although details must depend on global issues and boundary conditions so that one has 
to proceed with caution, we may also expect that this analysis gives reasonable order of 
magnitude estimates even when \kr\ ~ 1. 

In addition, for simplicity we shall also continue to consider the polytropic model as the 
essential physics is contained therein, while bearing mind that results, when expressed only 
in terms of the local disc model parameters, S, H, and p*Co; can be straightforwardly applied 
to other models on the basis of the discussion given in section 16.51 above. 

Setting X, W oc e*'^'" in equations ( !67l) and (!68l) we obtain 



respectively, where k = kr is the dimensionless wave number. One can use equations (17T1) 
and ( 1721) to obtain dispersion relation u = uj{k). However, as we have to consider several 
limiting cases for which there is qualitatively different dynamics, the general form has limited 
usefulness. 




(71) 



and 




(72) 
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7.1.1 Dynamics of a twisted disc without external radiation pressure: Viscous and wave 
regimes 

At first let us recall the basic dynamical properties of the standard twisted disc where radi- 
ation pressure effects are absent. In this case the dispersion relation can be easily obtained 
from equation (!72l) together with the condition x = 0. We then obtain 



u = -{-ia±^j6^kyQ-a^). (73) 

The behaviour of the disc is qualitatively different in the two limiting cases of sufficiently 
large and sufficiently small a parameter (PP, Papaloizou & Lin 1995). For the case of a large 
a the 'slow' mode corresponding to the the positive imaginary part in (1751) determines the 
disc dynamics and the dispersions relation can be approximately written as 

Co -iujy, = ——. (74) 

Remembering that according to our convention all quantities are proportional to e"""^* ~ 

^-uj„t conclude that in this case, twisted perturbations of the disc decay. In the limit of 

small a it is easy to see that twisted perturbations propagate over the disc with a speed of 

the order of a typical sound speed with little dissipation and we have 

6\~k\ , , 

cI; ^ ±uj, - ta/2, = ^-j=. (75) 

From the condition u^, < ojg one can see that the 'viscous' regime of evolution of the disc is 
realised when 

a > (76 
2^6 ^ ^ 

For the case of a 'global' perturbation with ~ 1 the above condition is approximately 
reduced to requirement that the viscosity parameter a is larger than the disc opening angle 
S. 

7.1.2 The effects of radiation pressure for relatively large a 

Now let us consider effects induced by external radiation pressure. At first let us assume that 
the condition (176!) is fulfilled. In this case we have from equation (!72l) that approximately 

x ^ (cD + (77) 



Substituting equation ( 1771) in ( 17Ti) we obtain 



1 / 16 
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where we now define the time scale characterising influence of the radiation pressure, uu^ ^, 
such that 



The dimensionless time scale u)"^ is that associated with the radiation pressure driven in- 
stability considered by Pringle 1996. 

It is clear from equation (!78|) that the character of the dynamical evolution in fact 
depends on a comparison of Ur to A = p^/pc- As Ur = e^^^, this leads directly to the 
condition <^ A being required for the assumption that the vertical response of the disc 
is that of a rigid body to be valid for disturbances with radial scale comparable to r. That 
condition was also derived on general grounds in section 16.51 and from simple arguments in 
the introduction and is explored further below. 

First let us assume the case uj,. ^ A. In this limit the factor in the front of the flrst 
brackets in ( 1781) is equal to one and because by necessity A < 1, the flrst term in the 
brackets can be neglected, and therefore a) is purely imaginary such that 



This gives either growth or decay of the perturbation depending on whether the sign of the 
expression in the brackets is positive or negative. Since Ur is proportional to k it can be, in 
principal, positive when A; > 0. Thus one can infer instability of the disc as has been done 
by Pringle 1996. The quantity A drops out and the disc behaves as a collection of rigid rings. 
The corresponding stability criterion can be obtained from the condition that the absolute 
value of the flrst term in the brackets is larger than the second term and therefore 

Stt \k\ 
64 a 

Note that the numerical factor, ~ 0.1, in the expression for ecru depends on the vertical 
structure of the disc. It seems reasonable to assume, however, that it is always smaller than 
one for any possible vertical distribution of pressure and density. 

Let let us now consider the limit ujj. ^ A. In this limit the character of the disc evolution 
changes drastically. This is essentially because the vertical motion can no longer be treated 
as uniform as would occur for a series of rigid rings. Instead of potential growth of disc 
perturbations, we flnd mainly decaying oscillatory dynamics. The corresponding dispersion 
relation reads 



(79) 




(80) 
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4 2 

= —\ + ii-—-Cu^). (82) 

OTT IT UJr 

As seen from fl82l) the real part of the frequency does not depend on the magnitude of 
the radiation pressure. The imaginary part can be either negative or positive. The values of 
the terms in the brackets, being always negative in the limit of very small A or very large Ur 
and thus very large Fq, lead to damping in those limits. In any event we might also expect 
that the oscillations in the disc inclination could result in strong self- shielding of the disc 
which can inhibit the growth of the inclination angle. 



7.1.3 The effects of radiation pressure for small a 

Let us now consider the opposite case of a small viscosity parameter a such that the in- 
equality ( 176|) is not satisfied. In this case, provided that the value of radiation flux is smaller 
than the typical mid plane pressure of the gas, a simple analysis shows that the effect of 
radiation pressure is to give a correction to the dispersion relation (1751) . 

We begin by pointing out that equation fl72]) in this limit can be brought to the form 

x^-cj2 + 1 W, 83) 

where the last term in the bracket is a small correction. After substitution of (1831) in equation 
fl7T|) we can look for solution for the dispersion relation in the form; uj = ±uJs ~ ia + uji. It 
is easy to see that the correction ui has the form 

It is clear from fl84l) than we always have < ^(^r- In addition we note that 

pr/cDsl ~ e6. (85) 

As follows from the definition of the parameter e (see equation flMl) ). the right hand 
side of f l85|) is smaller than unity provided that the radiation flux is smaller than the disc 
midplane pressure. Assuming that this is valid the disc dynamics, in the case of the small 
a, is mainly determined by propagation of waves with a speed of the order of the sound 
speed as in the case without radiation pressure. As in the case of very small A considered 
above, we again get damping in that limit. It also seems reasonable to suppose that in this 
situation, fast oscillations of the inclination angle would lead to a strong self-shielding of the 
disc thus inhibiting possible instabilities. 
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7.2 Physical character of the disc dynamics in the regime where A ^ cj^ 

As discussed above, when X <^ ujr the character of the disc dynamics differs drastically from 
that of a set of rigid rings. The disc then tends to oscillate at the frequency ci) ^ ^A, see 
equation (182!) and the discussion in section 16. 5[ Since this regime has not been discussed 
elsewhere we briefly describe its main features. 

Assuming the disc to be in this regime, as indicated in section 16.51 we can neglect unity 
in the brackets on the right hand side of equation flTTj) and obtain a simple relation between 
X and W such that 
4A 

X ^ — W. (86) 

Sir 

Now we use equation (!86|) and equation (!69|) together with equation (lAlip of Appendix |A] 
to see that the radiation pressure term F_(_ — > when X/ujr —>■ 0. 

Therefore, in this limit the disc tends to evolve in such a way that its surface becomes 
almost parallel to the initial or radial direction thus strongly decreasing the amount of any 
increase in the amount of intercepted radiation as a result of the perturbation. The disc 
gas does not cross this surface and therefore, in this limit, the radiation pressure effectively 
provides an impenetrable wall at the density p* thus prohibiting the gas motion in verti- 
cal direction. However, when there is no influence of the radiation pressure on the lower 
disc boundary, as mentioned already in section 16.51 stable radially propagating warp like 
disturbances can still still exist. 



7.3 Maximal potential growth/decay rate 

An interesting consequence of the effects described above is that, when equation (1751) is 
considered, there is an extremum of the imaginary part of uj, uj as a function of Ur- The 
extremum is attained when (D^ = a)^ = ±|A and the corresponding value of ui is equal to 

u'j = ±^X-u,. (87) 
on 

Thus for any value of the external radiation pressure and corresponding positive/negative 
uJr corresponding to the instability/decay of the disc perturbations, the growth/decay rate 
cannot exceed the value given by equation fl87|) with the signs (+) and (— ) corresponding 
to instability and decay, respectively. 



30 P. B. Ivanov and J. C. B. Papaloizou 

When the effect of the external radiation pressure is significant, the second term on 
the right hand side of ( 187|) can be neglected and we have a very simple expression for the 
maximal absolute value of the potential growth/decay rate 

Pll = ^A. (88) 

on 

7.4 Conditions on disc model parameters for potential instability 

As indicated above, the disc dynamics in the presence of the radiation field is more complex in 
the case of a relatively large a for which inequality (1761) is fulfilled. Then, different dynamical 
regimes are possible depending on the parameters of the unperturbed disc model as well as 
the magnitude of the external radiation flux. Here we estimate the characteristic parameter 
values that separate the different dynamical regimes assuming that inequality (176|) is valid. 

According a the linear theory of twisted disc perturbations that makes use of a local 
WKB analysis, instability associated with the effect of radiation pressure appears to be 
unaffected by the disc density p^: at the location it is applied when \LJr\ ^ |A. When this 
condition is not satisfied, the disc dynamics is strongly affected by effects associated with 
the non-rigid body character of the response of the disc vertical structure and the potential 
instability may be either inhibited or absent. 

The above inequality leads to an upper limit for the value e and, accordingly, on the 
value of the radiation momentum flux Fq 

^ < = Im 

On the other hand the quantity e should be larger than ecru given by equation (ISTl) . Therefore, 
instability seems to be possible only for a certain range of values of the radiation momentum 
flux. From the condition ecru < we find a condition on the disc parameters for the 
instability to be possible in the form 

Aa vr , , , , , , , , 

^>-|A;|^0.1|A;|. (90) 

Assuming that this estimate is approximately valid even when \k\ ~ 1 we obtain a lower 
bound for the value of the density ratio A from ( l90i) 



A = ^ > X^rit = ^, (91 

Pc "(-1) 

where 5(-2) = S/10~^ and tt(-i) = When the density ratio is very small and such 

that A < Xcrit, the simple theory indicates that the disc is likely to be stable for the whole 
range of possible values of the radiation flux Fq. 
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7.5 Some conditions for non-linearity 

We can obtain a value Fq = of the radiation momentum flux corresponding to the 
characteristic value separating two possible regimes of evolution of the disc in the presence 
of radiation pressure. By making use of equation (|6^ we find 



F, ^ yPe, (92) 


where Pc = pS^'^^Q estimates the gas pressure in the disc mid plane. As mentioned above, in 
a situation where fiaring of the disc is present, the magnitude of the radiation momentum 
fiux should be smaller than P^ in order for our simple theory to be approximately correct. 
This sets an upper bound for the value of A, such that \ < 5. Assuming that the radiation 
momentum fiux, Fq, is near to F*, the expression fl92|) leads to a stringent constraint on the 
value of the inclination angle (3. In order for the linear theory to be valid, a characteristic 
value of the radiation pressure in the disc F ~ FQr{d(3 / dr) should be smaller than the value of 
the gas pressure at the surface where the radiation pressure is applied, = P{p,) ^ PcX^/^ 
for the polytropic model with 7 = 5/3. 

From the condition F < P^, assuming that (3 changes on scale comparable to the radius, 
we obtain a condition for linear theory to be valid in the form 

p<p, = X^/^ < 1. (93) 

Note, however, that this constraint has a significant dependence on the vertical structure of 
the disc. Thus, for an isothermal disc where P (x p the small parameter A^^^ in (l93ll would 
be absent. 

It is also interesting to note that for a polytropic disc, similar considerations lead to 
a strong constraint on the value of (3 when the radiation fiux Fq is of the order of Fcru 
corresponding to onset of a radiation pressure induced instability in the disc. In this case we 
can repeat the above argument but using the value ecru given by equation (|8T]) with k = 
instead of e*. 

Proceeding this we find that linear theory would be valid only if 

/3 < /3,„t = (94) 

Substituting typical values 6 = 10^^ and a = 10^^ and assuming that A ~ 10^'^ we get Pcru = 
10^^. Thus, for unperturbed disc models with strong density and temperature drops toward 
the surface boundary of the disc, the linear theory of perturbations is, strictly speaking, 
valid only for quite small values of the inclination angle. 
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7.6 The one dimensional vertically stratified slab analogue 

As we discuss in more detail in the next section there is an analogy between the motion in 
this regime and that occurring in the simple hydrodynamic problem of free one-dimensional 
oscillations of a vertically stratified column of gas, with equation of state ([8]), gravitational 
acceleration g = —^I'^C, lower free boundary, and an upper rigid boundary at some density 
p* <^ Po- The unperturbed distributions of pressure and density are given by equation (|TOl) . 
Assuming that perturbed quantities depend on time through a factor oc e~*(^+'^)* the free 
oscillations of the gas column are described by equation (H5il and therefore expressions ( 155|1 
and fl56l) with the response Vg = provide the solutions of the free problem for a polytrope 
with 7 = 5/3. As above we shall focus on that case below. 

If we assume that such a column is rotating around a central body with angular frequency 
Q and m = 1, the frequency of oscillations observed in a non rotating frame will be equal to 
uj. As we will see later the smallest eigen value u, for the normal mode oscillations of this 
problem, is precisely equal to - the frequency of the disc oscillations, when the radial 
wavelength is very long, and it is in the regime X <^ Ur- 

8 LINEAR AND NON LINEAR DYNAMICS OF THE VERTICALLY 
STRATIFIED SLAB 

When the condition cD^ ^ A is fulfilled, internal non uniform vertical motions determine the 
disc dynamics. As we have mentioned in the previous section, the dynamics of the disc in 
the long radial wavelength limit is quite similar to the dynamics of a perturbed vertically 
stratified one-dimensional gas column with equation of state given by ([H]) and gravitational 
acceleration g = —Q'^(. 

Contrary to the disc case, the one dimensional problem can be easily studied, both 
analytically and numerically. The results can be used to check the validity of some of our 
assumptions leading to conclusions about the disc dynamics. They can also provide some 
understanding of the disc dynamics in the non-linear regime when the external pressure 
perturbation is of the order of disc pressure at the density level where the external forcing 
is applied. 

It turns out that the analogy between the disc problem and the one dimensional col- 
umn problem is practically exact in linear perturbation theory, provided that we consider 
perturbations of the column subject to the upper boundary condition 
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AP^P-Po = a{^)v, (95) 
"so 

where Pq is the equihbrium distribution of pressure given by equation (fTOj) . v is the gas 
velocity in the vertical direction, Pe is the externally applied pressure and a is a dimension- 
less parameter. But note that if unphysical behaviour is to be avoided with this boundary 
condition, we require a < 0, see below. All upper boundary quantities are assumed to be 
evaluated at some coordinate C* corresponding to the density p* = p(C = C*) ^ Pc- Note 
that in the limit a — oo, we must have f (C = C*) = in order to have a finite pressure at the 
boundary. Then the condition (|95|) reduces to a rigid wall condition. The lower boundary 
remains free. 

We go on to consider the dynamics of the vertically stratified column in detail, both 
analytically and numerically. As in the disc case we set 7 = 5/3. 

8.1 The eigenvalues for the linear modes of the one dimensional slab 

Let us assume that the perturbed quantities depend on time through a factor oc e~*°"*. In 
this case the free oscillations of the gas column are described by equation (H5l) for any value 
of a provided that we substitute the factor — for the factor 2Quj in the second term on 
the left hand side. Clearly, the solutions to the problem can be found from equations (153|) 



and (|55|1 with the parameter k expressed in terms of cr as K=yl + 3^. 



8.1.1 The case of two free boundaries 

If we assume that velocities are regular at both boundaries, taken to be where the density 
vanishes, this regularity condition determines a set of eigenvalues with frequencies cr„ ^ Q. 
The set of as well as the associated eigen functions can be easily determined from ([S5D- 
It is important for our discussion that the mode corresponding to the smallest eigenvalue, 
(T = has a velocity that is independent of Therefore this mode describes a uniform 
translation of the column as a whole. 



8.1.2 The case of an upper rigid boundary 



When we adopt the boundary condition fl95|) together with the regular boundary condition 
at C ~Co5 the eigen frequency of the shift mode acquires a small correction Q ^ Q + u. 



^ The regular modes with (T„ > H are just standing sound waves. 
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The resulting mode may be called as a 'modified shift' mode. As we mentioned above one 
can consider the column to be rotating around some centre with angular frequency Q. In the 
non-rotating frame, the shift mode of the free system describes a stationary displacement 
of the column as a whole, while for the case with a rigid boundary, the modified shift mode 
corresponds an oscillation of the column with characteristic time scale ~ l/cj. 

The change to the eigenfrequency, u, can be readily be determined by applying the 
boundary condition (jHSD to the solution given by equation flS^ under the assumption that 
\a\ and A = p*/pc ^ 1- 

This is found to be approximately given by 
H-^) 4a 



UJ 



5Aj 

;i + ^(f)^)15vr' 



(96) 



where we recall that uj = uj/VL. It follows from fl96l) that when a > the perturbations of the 
column decay with time. The condition a < leads to instability. Making the identification 
Cjr = —jqO. we see that equation fl96|) is precisely equivalent to equation f l78|) with Cjy = 0. 
Therefore, the linear dynamics of the gas column with boundary condition (!95|) should 
capture the form of the vertical motions relevant to the disc dynamics. 

We see from equation (l96ll that the real and imaginary parts of uj, namely ujr and uji are 
functions of only one variable y = a/X. Thus we have 
ujn ^ 4 uj_ ^ A y 

A TStt (1 + ^2/2)' X 15vr(l + ^|/2)- V ; 

Thus, all models with different a and A but the same ratio of these quantities have the same 
value of UJ expressed in units of A . As seen from equation fl971) the absolute value of uji has 



a maximum when y = h with \uji\ = j^. This corresponds to the condition (IHHIl . 



8.2 Non linear numerical solutions for the vertically stratified slab 

In our numerical work we use the simple Lagrangian staggered leap frog scheme described 
in Richtmyer & Morton 1967 for an ideal gas with 7 = 5/3. We use 400 grid points in the 
vertical direction uniformly distributed with respect to the vertical coordinate in the static 
initial state. Artificial viscosity is used in order to handle possible shocks and it is checked 
that our results are essentially independent on the value of the artificial viscosity coefficient. 
In order to check the scheme we compared the eigenvalues for the free problem (T„ and the 
Rankine-Hugoniot relations for shocks with what is obtained in numerical computations. In 
both cases we obtained good agreement between analytical and numerical approaches. 
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In our computations we use the dimensionless time r = Qt. The density p, pressure p, 
velocity v and position of a gas element are assumed to be expressed in units of pc, Pc, ^Co 
and Co; respectively. These normalisations are used to determine the normalisations of all 
other quantities that make them dimensionless. For example the thermal energy per unit 
area of the unperturbed disc, Efh = d(p. Following our procedure, this quantity is 

expressed in units of PcCo- Thus we have E^h = f| ~ 0.33 for 7 = 5/3 in these units. 

8.2.1 Numerical check of the eigenvalues and structure of the linear eigenmodes 

In order to check the validity of the linear theory we consider the case a > corresponding to 
decay of warp like perturbations. The case a < leads to a model with unphysical behaviour 
because it is such that, not only the modified shift mode, but also the higher order acoustic 
modes are unstable. The growth rates of the latter are always larger than that of the modified 
shift modes leading them to dominate the evolution. 

For a > 0, we start our computations with initial state having the density and pressure 
distributions corresponding to hydrostatic equilibrium and a very small uniform vertical ve- 
locity V = 10~^. We evolve the system for a very long time r ~ 10^ — 10^. Then we determine 
power spectra for the state variables and use them to locate the peak corresponding to the 
modified shift mode. It is confirmed that the associated spectral feature or line has a Lorentz 
profile. The location of the maximum determines the real part of the eigen frequency ur, 
and the width of the line determines uj. 

We consider four different values of A = 0.032, 0.056, 0.084, 0.11 and seven different values 
of a for a given value of A, such that a = 80X/2'', where the integer k = 1, 2. .7. When a is 
smaller than 5A - the value corresponding to the maximal value of ui we expect the column 
dynamics to be similar to the disc dynamics in the regime described by Pringle 1996. In this 
case the values of velocity and the Lagrangian displacement should not depend significantly 
on ( and the column response to the external pressure is similar to that of a rigid body. 
However, when a > 5A, we expect the behaviour to be similar to the disc dynamics in the 
regime described in section 17. 2[ In this regime the upper disc surface behaves as if it is in 
contact with rigid wall while the regions beneath participate in stable oscillations. In this 
case velocities and displacements should have a strong dependence on ( near the upper 
surface of the column. In a case of a very large a velocities and displacements should be 
close to zero at C ~ Co- 
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Figure 1. Numerically obtained values of the real and imaginary parts of the eigen frequency, lor and luj divided by the 
parameter A as functions of the ratio a/A. Solid curves represent analytical results given by equation l|97|l . Dashed and dotted 
curves represent numerical values of uiji and wi, respectively. The various symbol styles indicate different values of A. These 
are such that smaller values of A lie closer to the analytically determined curves. 

The results of our computations are plotted in Figures 1-4. In Figure 1 we show the 
numerically obtained dependencies of uji/X (the dashed curves ) and ujj/X (the dotted 
curves) on a/ A for the set of parameters described above. Specific symbols are associated with 
particular values of A. Circles, squares, diamonds and triangles represent the specific values 
of A we used in increasing order. Solid curves correspond to the analytically determined 
dependencies given by equation (197|) . One can see that the agreement between the analytic 
and numerical results is good for the smallest value of A = 0.032. However, the curves deviate 
more significantly as A is increased. This departure from a dependence of u/X only on a/A 
may be explained as an effect due to the increasing importance of of higher order terms in 
the power series representation of cj/ A as a function of A, for fixed a/ A, that were not taken 
into account in equations (!60l) and (!6T!) . 

Figure 2 shows the displacement vector evaluated near the boundary = (^o as a function 
of time for A = 0.032 with a = 1.26 (solid curve ) and a = 0.02 ( dotted curve). Regardless 
of the fact that these values of a differ, their analytically calculated decay rates are prac- 
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Figure 2. The time dependence of tlie displacement vector evaluated near the upper boundary of the gas column calculated 
for models with A = 0.032. The solid (dotted) curve corresponds to a = 1.26 (a = 0.02). Note that both curves practically 
coincide. The dashed curve represents the exponential dependence of the amplitude of the oscillations on time, being oc e~'^'*, 
calculated with the analytically determined value of u)/ = 1.65 • 10~^. 

tically the same and such that ui ^ 1.65 ■ 10^'^. As expected the temporal behaviour of 
the displacement vector is almost indistinguishable for these two cases. The dashed curve 
indicates the exponential dependence on time, calculated in this case, using the decay rate 
cui — 1.65 • 10~^. It is evident that agreement between the analytical and numerical results 
is good. 

Figure 3 is similar to Figure 2 but the dependence of the displacement vector on time 
is shown for a = 0.16.. Because a/A = 5 this value of a corresponds to the maximal decay 
rate when A = 0.032. In this case it is given by (x;/ ~ 6.7 • 10""^. Note again the very good 
agreement between the theoretical and numerical results. 

In Figure 4 we show the dependence of the displacement vector on ( calculated at the 
same time for a — 1.26 - solid curve, a — 0.16 - dot-dashed curve and a — 0.02 -dashed 
curve. One can see from this Figure that the displacement corresponding to the small value 
of a is practically independent of (. But for the large value of a it has a prominent increase 
toward the upper boundary of the column. Also, note that in the latter case the value 
of displacement vector near the upper boundary is close to zero. Thus, although the time 
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Figure 3. As Figure 2 but with a = 0.16. This value of a gives the maximal decay rate when A = 0.032. 



Table 1. The initial amplitude of velocity profiles, C, the maximal value of the relative pressure difference 
(P — Po)/Po = ^P/Po calculated for the first pulsational period and the difference between the final and 
initial values of the thermal energy per unit area, AE-p^, for the different models used in the computations. 
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dependence of the displacement vectors corresponding to the cases of large and small a is 
practically the same, the spacial dependence of these quantities is quite different. 

8.2.2 Non linear calculations 

In order to determine how our results depend on the initial velocity amplitude, we consider 
the case for which the upper boundary of the column is a rigid wall, resulting in the condition 
that f (C — C*) — 0- Prom our previous discussion it follows that this occurs when the 
parameter a is sufficiently large. The uniform initial velocity profile used in our previous 
numerical calculations is not compatible with this boundary condition. Therefore, we use a 
different initial form for the velocity 
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Figure 4. The dependence of the displacement vector on the vertical coordinate at a typical time, for A = 0.032 and different 
values of a. The solid, dot-dashed and dashed curves correspond to a = 1.26,0.16 and 0.02 respectively. 



V = Ccos{-—) 



(98) 



Since the non-linear effects are expected to operate on a short dynamical time scale of order 



of a few periods P = 2ttQ we follow the evolution for relatively small times 



end 



200. 



When the rigid wall boundary condition is applied and the linear theory is valid, the 
column dynamics is determined by the linear superposition of strictly periodic oscillations. 
Departures from this type of motion are determined by non-linear effects. One of the most 
important for our purposes is the possibility of shock formation near the upper boundary. 
When such shocks propagate through the column they irreversibly heat up the gas thus 
causing dissipation of the motion induced in the column. 

The strength of non-linear effects is obviously proportional to the induced amplitude C. 
In the first row of Table 1 we show the eight values of C used in our computations and the 



° Of course, shock formation in a realistic system would be quite different from what we consider in the text. Since the upper 
boundary of the disc is close to the surface of unit optical depth, in a realistic situation, the shocks would be radiative. However, 
it is plausible that our study indicates the range of disc parameters for which shocks may be expected. 




Figure 5. The dependence of density p on at four different times for model 7, which is characterised by the relative 
pressure amplitude C = 0.16. The solid, dashed, dotted and dash-dotted curves correspond to the times t = 2.5,5,7.5 and 10, 

respectively. 



corresponding models, enumerated by numbers from one to eight. In all cases we assume 
that the rigid boundary is situated at the value of C* corresponding to A = 0.032. 

A useful measure of non-linearity of the system is the ratio of the difference between 
the pressure at the upper boundary and its equilibrium value AP = P — Pq, to Pq. For 
definiteness we estimate the maximal values of this quantity during the first non-linear 
pulsation and display them in the second row of Table 1. As seen from this Table when 
AP/Pq < 1, this quantity has an approximate linear dependence on C, while for larger values 
this dependence is somewhat steeper. The third row represents the difference between the 
thermal energy (per unit area) at the end of the computations at Tend = 200 and its initial 
value Efh ~ 0.33, in units of PcCo- As expected this quantity is approximately proportional 
to the square of C. 
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Figure 6. The dependence of the entropy per unit area, S, for the six models (3 — 8). As the entropy increases while moving 
from the lowermost to uppermost curve, the models range from 3 — 8 with increasing values of the initial amplitude C. 



8.3 Shock development and the formation of a rcirefied hot atmosphere 

When AP/Pq < 1 the dynamics of the model is well described by the linear theory. In 
the opposite case, non-linear effects are substantial. The system exhibits a complicated 
pulsational behaviour and strong shocks are observed. This regime is illustrated in Figure 
5 where the results of calculations for model 7, with large initial amplitude C = 0.16, are 
presented. We show the density profiles at four different times r = 2.5 ( solid curve ), r = 5 
(dashed curve ), r = 7.5 ( dotted curve ) and r = 10 ( dot dashed curve ). The density 
distribution at r = 2.5 is quite similar to the equilibrium density distribution. When r = 5 
and T = 7.5 shocks are observed propagating from the upper boundary in the direction 
of negative (. These shocks heat the gas up, thus dissipating kinetic energy. Also note 
the density drop at the upper boundary when t — 5 and r = 10. We have checked that 
this density drop has a persistent character and the density at the boundary at r = Tend 
is significantly smaller than the initial density while the pressure does not exhibit such a 
significant change. 

Thus, the presence of shocks in the system leads to formation of a hot rarefied region 
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near the rigid boundary. It is not clear however, how this effect would operate in twisted 
discs. In the latter case the boundary is defined by the condition that the optical thickness 
is one. The position of this boundary may move when any density redistribution becomes 
significant. In our case with fixed upper boundary location, the presence of this hot low 
density region leads to the suppression of shock formation after a time corresponding to a 
few orbital periods, an effect that may not occur if the upper boundary moved in such a 
way as to maintain the relative pressure variations. 

The effect in our case is illustrated in Figure 6 where we show the dependence of entropy 
per unit area, S — J dC,p\og{p/ p^) as a function of time for the models with indices larger 
than 2, in the natural units PcCo- Curves attaining larger values correspond to larger ampli- 
tudes. One can see from this Figure that when the initial amplitude is sufficiently large, i.e. 
C ^ 0.04 and accordingly AP/Pq ^ 2.5, half of the entropy change is roughly determined for 
r < 10 — 15. (i.e. approximately during the first two orbital periods). Then shock production 
becomes less efficient due to the mechanism discussed above. 

8.4 Estimate of dissipation time scale 

As mentioned above, it is not clear to what extent the dynamics of the column in the strongly 
non linear regime is similar to the non-linear dynamics of twisted discs. One can assume, 
e.g. that in twisted discs, the boundary where external radiation pressure applies follows a 
certain density level. In this situation it seems reasonable to assume that the shock formation 
might persist at a rate similar to what is observed at initial times r < 10. In this case the 
disc perturbations could be damped on a time scale tni which can be roughly estimated from 
the results presented in Figure 1. 

For model 7, assuming the thermal energy production is of the order of AEth = 0.002, 
we can estimate the non linear decay time as t„; ~ 10{Eth/ AEth)fl~^ ~ 1.5 ■ lO^fi"^ where 
we recall that Eth ~ 0.3 in our units. Thus, under these assumptions the perturbations may 
decay in approximately 200 orbital periods. This rate is characteristically the internal energy 
content within one surface scale height per orbital period and becomes noticeable when the 
relative pressure variation AP/Pq near the upper boundary exceeds ~ 2.5. 

Let us stress, however, that although indicative, because of an obvious sensitive depen- 
dence on the density structure near the boundary, such an estimate cannot be rigorously 
justified within our model. Clearly more comprehensive future investigations are needed. 
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9 DISCUSSION AND CONCLUSIONS 

In this paper we have presented a self-consistent calculation of the response of a twisted 
disc to the action of a radiation pressure force acting in upper layers of the disc. This 
is assumed to be due to the interception and subsequent reradiation of radiation from a 
central source. The radiation pressure force is assumed to be applied at a single density level 
p=K corresponding to optical thickness unity. The degree of twist and warping is assumed to 
be small enough that linear theory can be used to calculate the response. 

The analysis of Pringle 1996 modelled the disc as a set of concentric rigid rings in a 
state of Keplerian rotation. These were assumed to communicate with each other through 
the exchange angular momentum because of the action of viscous forces. Up to now there 
has been no consideration of effects arising from the response of the disc vertical structure 
to the externally applied pressure. 

In this paper we have considered the warping and twisting of an accretion disc taking 
account the response of the disc vertical structure to an external radiation field assuming 
that the effect of self-shielding of radiation by the global twisted disc is not significant. We 
developed a description of the evolution of the disc in terms of a pair of equations governing 
the one dimensional evolution of the inclination as a function of radius and time (see section 
M and equations fHTl) . fl67|) and fl68l) ). This description extends earlier ones, so that the 
influence of radiation pressure on discs for which warps occur in the low viscosity bending 
wave regime, as well as the higher viscosity diffusive regime, may be considered. 

We found several qualitatively different dynamical regimes that may be realised in astro- 
physical discs. These are related to whether the character of the response of the disc vertical 
structure causes significant departures from what would be obtained for a rigid body. 

9.1 Conditions for a quasi-rigid response 

We found that to avoid such departures, the momentum flux due to radiation from the central 
source Fq, should be smaller than a critical value, F*, given by = {\/6)Pc, where 6 = C^o/r 
is the disc aspect ratio, Pc is the disc mid plane pressure and it has been assumed that the 
the ratio (A/(5) is small (see sections [6.51 and [Y.1.2p . Then, when Fq > Fcru ~ 0.l{S/a)Pc, 
warping instability of the disc becomes a possibility ( eg. Pringle, 1996). 

Thus for radiation driven warping of a disc that behaves like an ensemble of rigid rings, 
we flnd two requirements, namely that F^, > Fcru, together with Fq < F^,. These conditions 
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together imply that the ration of the surface to mid plane density should be sufficiently large 
and such that A > \crit ~ 0.1 5^ /a. 

9.2 Large external radiation momentum flux 

In the opposite limit of large Fq > max{F^, Fcru), consideration of the disc structure in the 
vertical direction is essential as the vertical displacement ceases to be uniform. The perturbed 
gas motion in the disc is mainly determined by the vertical component of velocity and the 
perturbed quantities tend to oscillate at a characteristic frequency u = (pcCo/2S)Ai7, with 
Q being the local Keplerian angular velocity. In this situation, vertical motion tends to be 
suppressed by the external pressure field on the irradiated side of the disc while warping 
motions persist in the bulk of the disc and on the shielded side of the disc (see sections 16.51 
andQ. 

The disc surface intercepting the radiation tends to align parallel to the radiation rays, 
thus decreasing the amount of intercepted radiation, while the inclination angles associated 
with the disc mid plane and the opposite free boundary of the disc oscillate, being approxi- 
mately equal to each other. In this limit the upper surface of optical thickness one plays the 
role of an impenetrable wall. Thus the development of instability of the inclination angle 
due to radiation pressure effects (e.g. Pringle 1996) would be inhibited in this limit. 

In principle, the presence of warping motions in this regime would be associated with 
variability on a large characteristic time scale T^h ~ X~^Tx, where and A are some 
'typical' values of the orbital period and the density ratio. 

9.3 Limits on the WKB growth rate 

Following Pringle 1996 we have considered possible instabilities using a WKB approach. 
As this neglects potentially important global effects and boundary conditions results are 
not definitive, nonetheless they should give a good indication of when instability could be 
possible. 

We find that when Fq < F^ the growth rate increases with Fq while in the formal limit 
-Fq ^ oo it approaches zero. Thus in the linear theory there is an upper limit for the growth 
rate of the instability of the inclination angle which is always ^ (pcCo/^)-^^ in the absence 
of viscosity (see section I7.3p . When present, viscosity acts to damp any instability at a 
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rate 0.16'^Q/a, leading again to the condition A > Xcrit ~ 0.16"^ /a for radiation warping 
instability to be feasible. 

9.4 Conditions for non-linearity 

Our results described above rely upon the applicability of linear perturbation theory. As 
discussed in section 17.51 the breakdown of linear theory is expected when the Lagrangian 
change of pressure induced in the surface layers becomes of the order of the unperturbed 
pressure even though the change of inclination angle could be very small. The corresponding 
constraints on the inclination angle are especially strong for disc models with a significant 
drop of temperature toward the boundaries of the disc. 

9.5 Non linear calculations for the one dimensional slab analogue 

A direct numerical approach to the issues discussed above is not yet feasible due to three- 
dimensional nature of the problem and the many physical processes involved. However, when 
vertical motions dominate, the problem becomes analogous to the one dimensional problem 
of vertical motions induced in a stratified gas column orbiting around a central source with 
Keplerian angular velocity. 

As discussed in Section |8] the one dimensional slab gives the same fundamental period 
of oscillation as obtained from the full disc theory when an appropriate boundary condition 
on the upper surface of the column is specified. The dependence of the eigenfrequency 
on theoretical parameters as well as the effective presence, in the limit of small A, of an 
impenetrable wall at the upper surface of the column have been checked numerically. Where 
they can be compared, agreement between our analytical and numerical results is very good. 

We also used the one dimensional model in order to investigate possible consequences of 
non-linear behaviour in the system. To do this we studied the motion of the slab ensuing 
from the imposition of a vertical velocity profile with varying amplitude. We found that when 
the ratio of the Lagrangian pressure perturbation to the local value of the pressure at the 
upper boundary of the disc, AP/P^, becomes larger than 1 — 10, strong shocks propagating 
downward into the column are observed (see section 1831) . In principle, these shocks may 
lead to non-linear dissipation of energy stored in the vertical motion and in section 18.41 we 
made a very crude estimate of a possible time scale of 200 orbits for AP/P^ ~ 10. However, 
this result must be checked in framework of a more sophisticated numerical approach which 
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takes into account physical processes neglected in this study, most notably the effects of 
radiation transfer in the upper layers of the column. 

9.6 Issues for future consideration 

In a fully self-consistent study the dynamical effects of radiation pressure must be studied 
together with effects determined by the radiation heating of the disc atmosphere. This can 
be done in the most convenient way within the framework of the one dimensional model 
discussed above. 

In a realistic disc model, where effects due to the flaring of the disc lead to the interception 
and scattering of radiation in the disc photosphere, when axisymmetric and unperturbed, 
radiation heating can significantly influence on or even determine the value of the density 
ratio A. 

This parameter, being the ratio of the density at the absorbing surface to the mid plane 
density defines the boundary between different dynamical regimes for a twisted disc. In this 
connection it is interesting to note that in certain vertical models of X-ray heated accretion 
discs, the ratio A can be quite small, being of order of 10~^ — 10~^, e.g. Jimenez-Garate et 
al 2002. In such models the new dynamical effects discussed in this Paper would certainly 
play an important role. 
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APPENDIX A: EXPRESSION FOR RADIATION PRESSURE ACTING ON 
THE DISC SURFACE 

Expressions for the radiation pressure acting on a surface element of a twisted disc (or, 
alternatively an expression for the corresponding torque acting on a disc annulus of radius 
r) have been derived by a number of authors, e.g. Petterson 1977b, Pringle 1996, Ogilvie & 
Dubus 2001, under different approximations. However, these authors assume that the disc 
consists of rigid rings of different radii having differing inclinations and orientations with 
respect to a fixed Gartesian coordinate system, when deriving the corresponding expressions. 
This assumption is essentially equivalent to the statement that all vertical motions can be 
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accounted for by a rigid tilt and it precludes internal vertical motions that vary with the 
vertical coordinate. 

As discussed in the main text, this is not strictly valid, as in general such motions occur. 
Accordingly, the density surfaces corresponding to an optical thicknesses of order of unity 
where the radiation pressure force is applied are not everywhere parallel to the surface locally 
corresponding to the disc mid-plane and the expression for the radiation pressure acting on 
the disc must be, accordingly modified. 

In general, an expression describing the surfaces of constant optical thickness can be 
easily obtained from the transformation law from the Cartesian coordinate system to the 
twisted coordinate system written for a gas element situated near the surface of the disc 

a; = r COS0 ± /3Co sin7, y = rsin0 =f /5CoCOS7, (Al) 

and 

z = ±Co + e^(C ^ Co) + r/3sin(0 - 7), (A2) 

where the upper (lower) sign corresponds to the upper (lower) disc surface. Let us recall 
that we assume that the inclination angle /3 is small and derive all equations in the linear 
approximation for the angle f3 and the Lagrangian displacement vector C,- 

The unit vector perpendicular to the disc surface is given by the standard expression 

(A3) 



_ _dR OR 
dr d(h 



dR OR 

dr dd) 



where R is the radius vector with components (x, z). 

We make the standard assumption that the radiation propagates radially before being 
intercepted by the disc and neglect the effect of self-shielding. In this case an amount of 
radiation intercepted by the disc and subsequently reradiated away is proportional to cosine 
of angle between the radius vector and the normal to the disc, cosS = n ■ R/\R\. For the 
upper branch of the surface a simple calculation gives 

d (Co + .... 

cos^ = —r— rW, (A4) 

or r 

where we recall that W = /5' sin-ip — P'y' cos ip = sin0 — \&2 cos0. As we discuss in the 
main text, in this paper, for simplicity, we neglect the effect of disc flaring and accordingly 
set |:(^) = in equation flAil) . 

The upper (lower ) branch of the disc surface is able to intercept the radiation of a central 
source only when cosH < ( cosS > 0). Therefore the radiation pressure, F, applied to the 
disc surface has the form 
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F± = T0(Tcos(S))^Focos(S), (A5) 

where B(x) is the step function, the signs (+) and (— ) stand for the upper and lower branches 
of the disc surface, respectively, 

" ^ (A6) 



47r cr^ 

is the momentum flux per unit of surface carried by the radiation emitted by the central 
source, and L is the luminosity of the source. The factor (2/3) accounts for (assumed ) 
isotropic character of the disc reradiation. 

In order to have the same symmetries as the surface pressure term (lASp the (^-component 
of the displacement vector, C,'^, and, accordingly, the (^-component of velocity perturbation, 
f^, must satisfy 

e'^(C,0) = -eH-C,0 + vr), v'^iCA) = -vH-CA + ^)- (A7) 

On the other hand, for a given value of (j), may be separated on even and odd components: 
f = + v'l and both components should be taken into account in order to satisfy the 
appropriate boundary conditions at C ^ ±^o- However, the latter requirement poses a 
problem with the symmetry law (]A7p . Indeed, while the even part of ^ f^, agrees with 
flA7p the odd one is not and the standard separation of the velocity perturbations into 
even and odd parts although being respected by the perturbed equations of motion, is not, 
strictly speaking, compatible with the symmetries provided by the boundary conditions. This 
problem can be circumvented with the help of the assumption that the normal decomposition 
into even and odd parts in C, of the perturbed vertical velocity is replaced by 

= — sign{sm{(l) — (A8) 

where ff and have even (odd ) properties with respect to the reflection ( Since 
the surface radiation pressure vanishes when cosH = or equivalently sin(0 — (pi) = and, 
accordingly, we have f ^ = at the corresponding values of = 0i and = 0i + vr the 
velocity anzatz ( lASp is self-consistent. Also, it agrees with the symmetries determined by 
equation (]A5I) and it has local simmetries respected by the equations of motion. 

Above and in the main text we always consider the angle to be in a segment [0i..0i + 
vr] defined by the condition sin(0 — 0i) > 0. With help of this convention we can treat 
our dynamical variables as having the standard even and odd symmetries with respect to 
reflection ( —(. Continuation of these variables into the segment defined by condition 
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sin(0 — 0i) < can be made with the help of the expression (]A8|) . In the case when 
sin(0 — (pi) > it follows from equations ( 1A4I) and ( lASj) that we have 

F{p = p.X^ Co) - = ^For|: + , (A9) 

and 

Fip = p,,C^ -Co) =F^ = 0, (AlO) 

where we recall that the density corresponds to the surface of an unperturbed disc with 
optical thickness r ^ 1. Taking into account that we can write = F| cos + sin and 
= 'Ci cos + ^2 sin 0, we can represent equation flAQI) in a convenient complex form 

where F+ = F| + iFl, = + and W = /Se*'^ = + z\l/2. Obviously, we also have 
F_ = 0. 



APPENDIX B: SOLUTION OF EQUATION (45) FOR THE CASE OF 
ARBITRARY 7 

When the parameter 7 > 1 is not specified, the homogeneous equation fH5l) can be reduced 
to a simpler form with help of equations ( fTOi) and ( fTTl) . We get 



2, d"^ 27X d ^ 4 



(1 - x^)— V - ^;^v + -^v = 0, (Bl) 

ctx^ 7 — 1 ctx 7 — 1 

where x = ^/^o- A general solution of (IBip has the form 

V = (1 - {CiF{^ -a,^-(3,^,x') + C2XF{1 - «, 1 - /?, | x')^ , (B2) 

where F{a, (3,S,x) is the hyper-geometric function and 



7 + 1 

a,p= J A \± 



1 + ^^). (B3) 



4(7-1)' ^ (7 + 1) 

. When ^ 1 we have 

7 + 1 8 8 

a^-^ + ^c:., (B4) 
2(7-1) 7 + 1 7 + 1 

APPENDIX C: THE GAUGE TRANSFORMATIONS 

As we pointed out in the main text different twisted coordinate systems can describe the 
same physical situation provided that they are connected by a certain family of transfor- 
mations (gauge transformations) which maintain fixed distributions of physical quantities 
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characterising the system (density and velocity fields in the case of a baratropic gas) in a 
non rotating (say, cylindrical) coordinate system with origin located at the central source. 
These transformation laws can be found by considering transformations from the cylindrical 
to the twisted coordinates. We derive them here assuming, for simplicity, that the equation 
of state is baratropic and the gas is inviscid. The latter assumption allows us to neglect the 
presence of the radial component of velocity in an unperturbed background state. 

We begin by considering the transformation of density. Let the density p be represented 
in the cylindrical and twisted coordinate systems as 

P = Poire, Zc) + Pi(rc, Zc) = po{r, Q + pi(r, C), (CI) 

where the background part pQ has the same functional dependence on the cylindrical and 
twisted coordinates. The cylindrical and twisted coordinates are related to each other by 
transformations ( fT6l) -( fT8l) . Substituting these transformations in equation ( ICip we obtain 

pUr,C) = Pi{r,C) + DPsmi^-j), (C2) 

where 

By construction the quantity pi is not changed after a gauge transformation. Therefore, it 
may be classified as a 'gauge independent' quantity. The result flC2l) can be represented in 
a compact form with help of the complex notation introduced in equation ( fT9l) 

p'i = Pi+ iDW. (C4) 

In an absolutely analogous way we can find the law of transformation of velocity pertur- 
bations. We have 

v<^- = v'^ + C(i(rfi)'W- W), (C5) 
= + irW + rfiW, (C6) 

and 

v*^^ = - <W - C^^W. (C7) 

Here the quantities v'^'^, and represent the velocity perturbations in the cylindrical 
coordinate system. They are gauge independent by construction. The quantities v"^, and 
v'' represent the time derivative of the corresponding coordinates of a comoving gas element. 
As we discussed above these quantities do not coincide with the projections of velocity onto 
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the orthonormal basis used in the formahsm exploiting the twisted coordinate system. The 
relations between these projections and the time derivatives of the coordinates can be easily 
found from the results provided in the Appendix of Petterson 1978, his equations (A3). They 
can be written in the form 

v-^ = v"^ - CW, = v"^ + irW, V = V - iCW, (C8) 

where it is implied that the background parts of the {() and (r) -components of velocity can 
be neglected. Substituting equations (ICSp in equations (IC5P -f lC7|) we find 

^<Pc = + <(rn)'W, (C9) 
v"^ = + rfiW, (CIO) 
and 

Now let us suppose that there are two twisted coordinate systems with two different Wi 
and W2 describing the same physical system. That means that the quantities characterising 
the density perturbation and perturbations of the velocity field in the fixed cylindrical system 
must be the same when expressed in terms of Euler angles and corresponding perturbations 
belonging to the different twisted coordinate systems. From this requirement and equations 
(lC4p and (lC9P - (]Clip it follows that the differences in the corresponding to these systems 
perturbations of density and velocity components must be proportional to the difference 
between Wi and W2, 5W = Wi - W2, 

5pi = -tD6W, (C12) 
= -iC{rQ)'SW, (C13) 
5v« = -rfi5W, (C14) 
and 

^v" = C^^^W. (C15) 
Equations (IC12l) - flC15P define the gauge transformations between two twisted coordinate 



system for an inviscid disc. Viewed in these terms, the generalisation to the case of a viscid 
disc is straightforward provided that the background part of the (r)-component of velocity 
is taken into account. 

One can also check that equations (1C12|) - (]C15|) provide an exact solution of the equa- 
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tions of motion and the continuity equation for an inviscid disc being written in a twisted 
coordinate system, for an arbitrary dependence of W on t and r. This solution describes 
the so-called trivial modes corresponding to an unperturbed disc. In this case the density 
and velocity fields in the cylindrical coordinate system remain equal to their unperturbed 
values while the dynamics in the twisted coordinates is such that perturbations of density 
and velocity defined with respect to the twisted coordinates are exactly compensated by 
those associated with time dependent Euler angles. 

It is important to stress that our dynamical equations ( |23l) and (1251) as well as the 
continuity equation ( l29i) are not invariant under the transformations generated by ( 10121) - 
f lClSP even when the term determined by viscosity are discarded. This is due to the fact 
that certain terms are neglected in these equations. These terms provide only next order 
corrections in our expansion series in small parameters when the gauge fl57j) is fixed. On 
the other hand the combination C.'^/r + zW entering in the expression (1A11I) determine 
an invariant quantity - the angle between the normal to the disc surface and the radial 
direction. This combination must be the same in all twisted coordinate systems connected 
by transformations (]C12|) - (]C15p . Using equations (1371) and (ICSP it is easy to see that this 
combination is proportional to 

v^ + rfiW. (C16) 

From equation (10140 it follows that the expression (]C16|) is, indeed, a gauge invariant quan- 
tity. 
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